Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
Revisit the analysis of He et al. (2010).
To get started, we must install pomp. Run the following to do so.
require(devtools)
install_github("kingaa/pomp")
set.seed(594709947L)
require(ggplot2)
require(plyr)
require(reshape2)
require(magrittr)
require(pomp)
stopifnot(packageVersion("pomp")>="0.70-1")
load("twentycities.rda")
measles %>%
subset(town=="London" & year>=1950 & year<1964) %>%
mutate(time=(julian(date,origin=as.Date("1950-01-01")))/365.25+1950) %>%
subset(time>1950 & time<1964, select=c(time,cases)) -> dat
demog %>% subset(town=="London",select=-town) -> demog
dat %>% ggplot(aes(x=time,y=cases))+geom_line()
demog %>% melt(id="year") %>%
ggplot(aes(x=year,y=value))+geom_point()+
facet_wrap(~variable,ncol=1,scales="free_y")
demog %>%
summarize(
time=seq(from=min(year),to=max(year),by=1/12),
pop=predict(smooth.spline(x=year,y=pop),x=time)$y,
birthrate=predict(smooth.spline(x=year+0.5,y=births),x=time-4)$y
) -> covar
plot(pop~time,data=covar,type='l')
points(pop~year,data=demog)
plot(birthrate~time,data=covar,type='l')
points(births~year,data=demog)
plot(birthrate~I(time-4),data=covar,type='l')
points(births~I(year+0.5),data=demog)
rproc <- Csnippet("
double beta, br, seas, foi, dw, births;
double rate[6], trans[6];
// cohort effect
if (fabs(t-floor(t)-251.0/365.0) < 0.5*dt)
br = cohort*birthrate/dt + (1-cohort)*birthrate;
else
br = (1.0-cohort)*birthrate;
// term-time seasonality
t = (t-floor(t))*365.25;
if ((t>=7&&t<=100) || (t>=115&&t<=199) || (t>=252&&t<=300) || (t>=308&&t<=356))
seas = 1.0+amplitude*0.2411/0.7589;
else
seas = 1.0-amplitude;
// transmission rate
beta = R0*seas*(1.0-exp(-(gamma+mu)*dt))/dt;
// force of infection
foi = beta*pow(I+iota,alpha)/pop;
// white noise (extrademographic stochasticity)
dw = rgammawn(sigmaSE,dt);
rate[0] = dw*foi/dt; //infection rate (stochastic)
rate[1] = mu; // natural S death
rate[2] = sigma; // rate of ending of latent stage
rate[3] = mu; // natural E death
rate[4] = gamma; // recovery
rate[5] = mu; // natural I death
// Poisson births
births = rpois(br*dt);
// transitions between classes
reulermultinom(2,S,&rate[0],dt,&trans[0]);
reulermultinom(2,E,&rate[2],dt,&trans[2]);
reulermultinom(2,I,&rate[4],dt,&trans[4]);
S += births - trans[0] - trans[1];
E += trans[0] - trans[2] - trans[3];
I += trans[2] - trans[4] - trans[5];
R = pop - S - E - I;
W += (dw - dt)/sigmaSE; // standardized i.i.d. white noise
C += trans[4]; // true incidence
")
dmeas <- Csnippet("
double mn = rho*C;
double v = mn*(1.0-rho+tau*tau*mn);
double tol = 1.0e-18;
if (cases > 0.0) {
lik = pnorm(cases+0.5,mn,sqrt(v)+tol,1,0)-pnorm(cases-0.5,mn,sqrt(v)+tol,1,0)+tol;
} else {
lik = pnorm(cases+0.5,mn,sqrt(v)+tol,1,0)+tol;
}
")
rmeas <- Csnippet("
double mn = rho*C;
double v = mn*(1.0-rho+tau*tau*mn);
double tol = 1.0e-18;
cases = rnorm(mn,sqrt(v)+tol);
if (cases > 0.0) {
cases = nearbyint(cases);
} else {
cases = 0.0;
}
")
initlz <- Csnippet("
double m = pop/(S_0+E_0+I_0+R_0);
S = nearbyint(m*S_0);
E = nearbyint(m*E_0);
I = nearbyint(m*I_0);
R = nearbyint(m*R_0);
W = 0;
C = 0;
")
toEst <- Csnippet("
Tmu = log(mu);
Tsigma = log(sigma);
Tgamma = log(gamma);
Talpha = log(alpha);
Tiota = log(iota);
Trho = logit(rho);
Tcohort = logit(cohort);
Tamplitude = logit(amplitude);
TsigmaSE = log(sigmaSE);
Ttau = log(tau);
TR0 = log(R0);
to_log_barycentric (&TS_0, &S_0, 4);
")
fromEst <- Csnippet("
Tmu = exp(mu);
Tsigma = exp(sigma);
Tgamma = exp(gamma);
Talpha = exp(alpha);
Tiota = exp(iota);
Trho = expit(rho);
Tcohort = expit(cohort);
Tamplitude = expit(amplitude);
TsigmaSE = exp(sigmaSE);
Ttau = exp(tau);
TR0 = exp(R0);
from_log_barycentric (&TS_0, &S_0, 4);
")
dat %>%
pomp(t0=with(dat,2*time[1]-time[2]),
time="time",
rprocess=euler.sim(rproc,delta.t=1/365.25),
dmeasure=dmeas,
rmeasure=rmeas,
toEstimationScale=toEst,
fromEstimationScale=fromEst,
initializer=initlz,
covar=covar,
tcovar="time",
zeronames=c("C","W"),
statenames=c("S","E","I","R","C","W"),
paramnames=c("R0","mu","sigma","gamma","alpha","iota",
"rho","sigmaSE","tau","cohort","amplitude",
"S_0","E_0","I_0","R_0")
) -> m1
m1 %>% as.data.frame() %>%
melt(id="time") %>%
ggplot(aes(x=time,y=value))+
geom_line()+
facet_grid(variable~.,scales="free_y")
readRDS("He2010_mles.rds") %>%
subset(town=="London") -> mle
paramnames <- c("R0","mu","sigma","gamma","alpha","iota",
"rho","sigmaSE","tau","cohort","amplitude",
"S_0","E_0","I_0","R_0")
coef(m1) <- unlist(mle[paramnames])
require(foreach)
require(doMC)
registerDoMC()
foreach(i=1:4,
.packages="pomp",
.options.multicore=mcopts) %dopar% {
pfilter(m1,Np=10000)
} -> pfs
logmeanexp(sapply(pfs,logLik),se=TRUE)
## se
## -3802.933702 0.492356
m1 %>%
simulate(nsim=10,as.data.frame=TRUE,include.data=TRUE) %>%
ggplot(aes(x=time,y=cases,group=sim,color=(sim=="data")))+
guides(color=FALSE)+
geom_line()+facet_wrap(~sim,ncol=2)
m1 %>%
simulate(nsim=100,as.data.frame=TRUE,include.data=TRUE) %>%
subset(select=c(time,sim,cases)) %>%
mutate(data=sim=="data") %>%
ddply(~time+data,summarize,
p=c(0.05,0.5,0.95),q=quantile(cases,prob=p,names=FALSE)) %>%
mutate(p=mapvalues(p,from=c(0.05,0.5,0.95),to=c("lo","med","hi")),
data=mapvalues(data,from=c(T,F),to=c("data","simulation"))) %>%
dcast(time+data~p,value.var='q') %>%
ggplot(aes(x=time,y=med,color=data,ymin=lo,ymax=hi))+
geom_line()+geom_ribbon(alpha=0.2)
bake <- function (file, expr) {
if (file.exists(file)) {
readRDS(file)
} else {
val <- eval(expr)
saveRDS(val,file=file)
val
}
}
require(foreach)
require(doMC)
registerDoMC()
bake("mif1.rds",{
foreach(i=1:15,.inorder=FALSE,
.packages=c("pomp","magrittr"),
.combine=c,
.options.multicore=mcopts) %dopar% {
mif2(m1,Nmif=50,Np=1000,
cooling.fraction.50=0.2,
transform=TRUE,
rw.sd=rw.sd(R0=0.01,iota=0.01,
rho=0.01,sigmaSE=0.01,
tau=0.01,amplitude=0.01,
S_0=ivp(0.02),E_0=ivp(0.02),
I_0=ivp(0.02),R_0=ivp(0.02)))
}
}) -> mfs
matplot(sapply(mfs,conv.rec,"loglik"),type='l',ylim=c(-3900,-3800),
ylab="loglik",xlab="mif2 iteration")
foreach(mf=mfs,.inorder=TRUE,
.packages="pomp",
.options.multicore=mcopts) %dopar% {
pfilter(mf)
} -> pfs2
print(sapply(pfs2,logLik))
## [1] -3921.309 -3803.307 -3835.341 -4758.305 -3809.474 -3816.246 -3820.660
## [8] -3817.373 -3811.493 -3819.388 -3818.031 -3859.014 -3829.719 -3816.378
## [15] -3810.664
He, D., E. L. Ionides, and A. A. King. 2010. Plug-and-play inference for disease dynamics: Measles in large and small populations as a case study. Journal of The Royal Society Interface 7:271–283.